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Abstract 

Nonlinear panel methods have no proof for the 
existence and uniqueness of their solutions. The 
convergence characteristics of an iterative, non- 
linear vortex- lattice method are, therefore, care- 
fully investigated. The effects of several param- 
eters, including 1) the surface-paneling method, 

2) the integration method of the trajectories of 
the wake vortices, 3) vortex-grid refinement, and 
4) the initial conditions for the first iteration 
on the computed aerodynamic coefficients and on the 
flow field details are presented. The convergence 
of the iterative-solution procedure is usually 
rapid. The solution converges with grid refinement 
to a constant value, but the final value is not 
unique and varies with the wing surface-paneling and 
wake-discretization methods within some range in the 
vicinity of the experimental result. 


Nomenclature 
AR = aspect ratio 

a ki '^kj A" influence coefficients defined in Eq. (4) 

c k J 

a,b ■ semi-axes of the ellipsoid in Fig. 3 

b = wing span 

= span bending-moment coefficient 
C Dl = induced-drag coefficient = D^/qS 

= lift coefficient = L/qS 
C = section lift coefficient 

C M = pitching-moment coefficient = M/qSc^ 

C N0R = norma l - f° rce coefficient = NOR/qS 

0^ = pressure coefficient = (p - P ra )/q 

ACp = local load coefficient on the wing 

c = 1 airfoil-section local chord 

c ma = mean aerodynamic chord 
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c r «* root chord 

D^ = induced drag 

F = vector of singularity intensities 

L = lift 

t = length of rolled-up wake 

M = pitching moment 

NOR = normal force 

N g ■ number of source panels on the body 

N^ = number of vortex panels on the wings 

n = unit vector normal to the surface, point- 

ing into the flow 

•i 

n^ = spanwise number of vortex panels 

n c = chordwise number of vortex panels 

p = static pressure 

q ■ free-stream dynamic pressure 

q^ = strength of source panel number (i) 

R = underrelaxation coefficient (Eq. (8)) 

S = wing area (wetted surface) 

u,v,w ■= perturbation velocity components 

V ro = free-stream velocity vector 

V = vector of the free-stream velocity com- 

ponents that are normal to the panels 

x,y,z = Cartesian coordinates in the downstream, 
right-wing, and upward directions, 
respectively. 

Ax ,Ax .= length scales in the wake and on the 
P wing panels, respectively 

a = angle of attack 

Tj = strength of bound vortex number (j) 

U = free-vortex shedding angle & 

<J> - = perturbation-velocity potential 

u “weighting parameter in second-order 

integration (Eq. (8)) 
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Introduction 

' Vortex-lattice methods (VLM) have been used to 
compute the lift of. wings for the last 35 years. 

They were an extension of Prandtl's lifting-line 
theory to lifting surfaces. Falkner'a works, 1-3 
probably the earliest in this area, laid the foun- 
dations for the method, using distributed vorticity. 
Discretization into a vortex lattice was introduced 
by Young and Harper, 1 * and many computer programs, 
solving for the unknown vortex intensities, were 
developed during the following years (e.g., Ref. 5) 
and applied to various problems. 

The VLM is basically one of many panel methods 
that are more accurately described as surface singu- 
larity methods. Generally speaking, panel methods 
distribute panels of singularities (sources, 
doublets, or vortices) on the surfaces of bodies 
that are immersed; in a flow. The strengths of the 
singular elements are determined by satisfying the 
tangency boundary condition on the body surfaces, 
and once these are known the influence of the body 
on any point in the flow field can be computed. 

Hunt 6 shows that the VLM is equivalent to a 
piecewise-constant doublet distribution and classi- 
fies it as a zero-order panel method. 

As such, and as long as compatible boundary 
conditions are used, 6 a unique solution exists for 
the vortex distribution. This can be proven, how- 
ever, only' for the so-called "linear" case, when the 
geometry of the wake that is generated by the vor- 
ticity distribution is known ab initio. 6 The classic 
example of a linear case is the streamwise-oriented 
planar wake of the finite wing in Prandtl's lifting- 
line theory. In practical cases the real wake is 
distorted after being shed and rolls up about its 
streamwise edges. Its geometry is part of the 
unknown flow field. The "nonlinear" panel methods 
are composed of complicated iterative techniques 
that "relax" the wake to its rolled-up, force-free 
shape and account for the two-way interaction of the 
lifting surfaces with their wakes by correcting the 
singularity distribution with every new shape of 
the wake. 

Nonlinear VLMs were developed in the last 
decade and were used extensively to compute the vor- 
tex lift of slender wings at high angles of 
attack. 7-18 These methods were extended to lifting 
slender bodies 19 ’ 20 and to high-angle-of-attack 
aerodynamics of multiple-wing body configurations of 
missiles and fighter aircraft. 21 Each step of the 
nonlinear iterative process in these methods is, in 
itself,' a linear solution. It recomputes the inten- 
sities of the vortices for every corrected wake 
shape using a prescribed wake geometry and has, 
therefore, a unique solution. However, there is no 
proof of the existence of a unique solution for the 
nonlinear procedure. This means that even when the 
iterative solution converges, one has no assurance 
of its uniqueness. This lack of proof of unique- 
ness motivated the present authors to conduct a 
careful study of the convergence characteristics of 
the method and its sensitivity to various geometri- 
cal and numerical parameters during the development 
of their nonlinear vortex-lattice computer program 
for wing-body configurations. 21 This paper describes 
the method briefly and presents its convergence 
characteristics vis-a-vis various parameters such as 
the vortex-panel size and shape, or the length of 
the free-vortex discrete segments. 


The Mathematical Model 

The method described in Ref. 21 was developed 
to compute the flow fields over winged flight 
vehicles at high angles of attack, and to predict 
their aerodynamic coefficients. The results of 
these computations for typical missile and fighter- 
aircraft configurations are presented there, and the 
full details of the numerical schemes are given in 
Ref. 22. Following is a brief description of the 
mathematical model used in this method. 

Assuming a steady, incompressible, inviscid, 
and irrotational flow (except for discrete vor- 
tices), the disturbance flow field of a vehicle 
flying in a uniform flow can be described by the 
Laplace equation: 

V 2 <f> - 0 (1) 

If the flow is subsonic and small perturbations are 
assumed, the following discussion will also apply 
when the Prandtl-Glauert equation is substituted 
for Eq. (1). 

The solution to the Laplace equation is deter- 
mined by the tangency boundary condition 

3<t>/9n + V„ . S = 0 (2) 

that has to be satisfied everywhere on the config- 
uration surfaces, and by the requirement that the 
perturbations vanish at infinite distances from the 
configuration. In the case of lifting configura- 
tions, their wakes are considered part of the con- 
figuration surface, 6 and the tangency boundary 
condition (Eq. (2)) has to be satisfied on them as 
well. This means that no pressure or velocity jump 
can exist across the wake or, in other words, that 
the wake must follow the local streamlines. This 
unknown location of part of the computational- 
domain boundary, which has to be determined itera- 
tively, makes the problem nonlinear, although the 
Laplace equation is linear. 

When Green's third identity is used, the prob- 
lem of finding the volume distribution of a poten- 
tial function is replaced by the problem of finding 
a surface distribution of singular elements 
(sources or doublets) that satisfies the boundary 
conditions. In the present method, the body, which 
would be nonlifting in potential flow, is modeled 
by a distribution of source panels only. The lift- 
ing surfaces are assumed to be of zero thickness 
and are modeled by a vortex lattice, where the dis- 
crete vortices are equivalent to a piecewise con- 
stant doublet distribution. Discrete free line 
vortices are shed from all the free edges of the 
lifting surfaces, and their trajectories in the 
wake are determined by integration of the local 
streamline equations 

dy/dx = v/(V TO + u) ; dz/dx = w/ (V„, 4- u) (3) 

from the shedding point into the far wake. At some 
predetermined distance downstream, the free vortex 
lines are replaced by semi-infinite straight lines 
in the direction of V„. 

Since any singularity distribution identically 
satisfies the Laplace equation (Eq. (1)) and the 
boundary condition at infinity, the simultaneous 
satisfaction of Eqs. (2) and (3) should generate the 
correct solution. Equation (2) must, however, be 
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satisfied on the wake also, the position of which, 
being part of the solution, is still unknown. 
Therefore, the solution procedure has to be itera- 
tive or have an impulsive start. In the iterative 
technique used here, the intensities of the singu- 
larities are determined by satisfying the boundary 
condition (Eq. (2)) on the surfaces of the config- 
uration and on an assumed wake shape. A new wake 
shape is determined by integrating Eq. (3), using 
the velocities induced by the previously computed 
singularities, and is used to recompute the intensi- 
ties of the singularities. This process is repeated 
until the solution converges. 


Method of Solution 

The body flow field is simulated by a conven- 
tional source-panel method. 23 The body surface is 
approximated by a number (N s ) of trapezoidal panels 
carrying a piecewise-constant source distribution 
of unknown strength q^ . The collocation point at 
which the boundary condition (Eq. (2)) is satisfied 
is at the center of the panel area. A typical 
paneling scheme of the body is shown in Fig. 1. 

The flow field generated by the lifting sur- 
faces is described by a nonlinear VLM. 18 The zero- 
thickness lifting surfaces are divided into a number 
(N v ) of vortex panels. These' are mostly rectangular. 
However, the triangular panels of Rom and Zorea 8 are 
used along the edges of swept: wings (Fig. 1). It 
will be shown later that this choice of the Rom- 
Zorea paneling had a first-order effect on the com- 
putational results. A discrete horseshoe vortex of 
unknown strength Tj is bound to the 1/4-chord 
line of each panel. The trailing arms of the vor- 
tex are also bound to the surface along the panel 
streamwise demarkation lines until they reach one 
of the wing edges. From there they are shed as free 
vortices (Fig. 1). The collocation points on the 
vortex panels are located at the intersection of the 
longitudinal median of the panel with its 3/4-chord 
line. 

A complete geometric definition of the problem 
also needs a specification of the wake shape. If 
better information is not available, the trailing 
vortices can be assumed to be semi-infinite straight 
lines leaving the wing edges at some predetermined 
angle above the wing surface; Past experience 21 has 
shown that any angle between the chordwise and the 
streamwise directions would do, but the choice of 
half the angle of attack is common. 

With the geometry defined, one has only to 
satisfy the tangency boundary condition (Eq. 2) 
simultaneously at all the collocation points. This 
results in a system of (N s + N v ) linear algebraic 
equations for the unknown intensities of the sources 
and vortices 

Ns N v 

2 a ki q i + 2 b kj r j + c k = 0 

1=1 j=i 

■' (k = 1, 2 N s + N y ) (4) 

where c^ = (V„ • n)j,. is the component of the free- 
stream velocity that is normal to the panel number 
(k) at its collocation point, and and bkj are 

the influence coefficients of the source panel (1) 
and the vortex (j), respectively, on the collocation 
point of panel (k). An influence coefficient is 


defined as the normal component (dQ/Bn)^ of the 
velocity that is induced on panel (k) by a singular 
potential element of unit strength. The influence 
of a vortex includes the influence of its trailing 
free arms, requiring the predetermination of their 
spatial location. Detailed formulae for the influ- 
ence coefficients are given in Ref. 22, and are not 
needed here for a discussion of the convergence 
characteristics . 

Rewritten in matrix form, Eq. (4) becomes 

[BW]F + V n = 0 (5) 

where [BW] is an (N s + Ny) square matrix of the 
influence coefficients, F is the vector of singu- 
larity intensities 



and V n is the vector V n = {c^}. The matrix [BW] 
can be partitioned into 



where BOB and BOW represent the influence coeffi- 
cients of the body sources on the collocation points 
of the body panels and wing panels, respectively, 
and WOB and WOW represent the influence coeffi- 
cients of the wing vortices (including the trailing 
free vortices) on the body panels and wing panels, 
respectively. BOB and BOW are computed only once 
because of the fixed geometry of the surface panels, 
but WOB and WOW have to be recomputed whenever the 
wake shape changes. 

Equation (7) is solved for the intensities of 
the singular elements (qi.Tj). This is a linear' 
problem and has, in principle, a unique solution. 
With these intensities known, the velocities induced 
by them on any point in the flow field can be com- 
puted and Eq. (3) can be integrated for the trajec- 
tories of the free vortices. The trajectories 
define a new wake shape that is used to recompute 
the singularity intensities from Eq. (7). The 1 : 
calculation cycle (Fig. 2) is reiterated until the 
wake solution converges to a constant shape (within 
a given tolerance) . 

Once the computation is converged, the load 
distribution on the lifting surfaces is computed 
from the Kutta-Joukowoski theorem, and the pressure 
distribution on the body is obtained from the 
Bernoulli equation. These distributions are inte- 
grated to compute the lift, induced-drag, and 
pitching-moment coefficients. Details of all the- 
calculations are given in Ref. 22. Only those 
necessary for the study of the convergence charac- 
teristics of the method will be discussed here. 

Validation of the Computer Code 

The objective of this investigation was to / 
study the convergence characteristics of the non- 
linear computation scheme for which there is no 
proof of the existence and uniqueness of the solu- 
tion. First, the computer code had to be validated 
in its linear mode to ensure its proper functioning 
when a unique solution existed by comparing the 
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results with other exact solutions or with proven 
numerical results. ' 

The source-panel body-simulation program was 
validated by comparing its results with the exact, 
analytically calculated, zero-incidence pressure 
distribution over an ellipsoid. Figure 3 depicts 
the relative error (AC p /C p ) in the pressure coeffi- 
cient as calculated by the present method. A mesh 
refinement (to 240 panels) resulted in a highly 
accurate solution except for a maximum error of 4% 
in the stagnation region, where a higher panel 
density would be needed to account for the large 
local gradients. . Accurate results (Fig. 4) were 
obtained on an ogive-cylinder at a .*■ 6° because 
of the smaller gradients. They are validated by 
comparison with the commonly used Woodward method. 24 

The vortex-lattice wing-simulation part of the 
program was validated first by comparison with 
analytical results for the chordwise circulation 
distribution on a two-dimensional flat-plate air- •' 
foil (at a = 4°) and for the spanwise circulation 
distribution oyer Prandtl's lifting-line represen- 
tation of a finite wing. The agreement in both 
cases was good. Figure 5 shows the relative error 
in the results for the flat-plate airfoil and the 
effects. of mesh refinement on these results (the 
results for the lifting line are not shown here). 

The computed circulation is within 0.25% of the 
theoretical value over most of the chord, even with 
only 10 vortex panels. The error increases near the 
leading edge because of its strong singularity and 
near the trailing edge where the denominator in the 
relative error goes to zero. It is interesting to 
note that the influence of mesh refinement (above 
N v = 10) on the integrals of the circulation (the 
section lift and pitching-moment coefficients) is 
insignificant 22 (not shown). 

Further validation of the VLM code was obtained 
on the delta wing, depicted in Fig. 6, with and 
without the presence of the body. Chordwise and 
spanwise load distributions at several locations on 
the wing, obtained with the Rom-Zorea paneling 8 
(Fig. 6a), were compared 22 with the results of the 
classic VIM trapezoidal paneling 5 (Fig. 6b) and 
with results of a trapezoidal paneling (Fig. 6c) 
with a piecewise constant doublet distribution. 
(Jacobs, A., Private communication, IAIFLO computer 
code, Courtesy of the Israel Aircraft Industry, 
Ben-Gurion Int'l. Airport, Lod, Israel, 1982.) One 
representative chordwise load distribution on the 
wing at the 20.5% span station (a = 6°) is shown in 
Fig. 7. These data were obtained in the presence of - 
the body that was modeled in all three programs by 
the same geometrical source-panel distribution. 

The continuity of the vortlcity distribution at the 
wing-body intersection line was obtained by the con- 
tinuation of the wing-root panels through the body 22 
(Fig. 6). The good agreement between the three 
methods in Fig. 7 and in Ref. 1 22 validated the 
linear mode of the present program. 


Parameters Influencing Convergence 

The solution to the nonlinear problem is 
obtained by an iterative procedure that is termi- 
nated when convergence is reached. One of the 
following three convergence criteria can be chosen 
arbitrarily according to the user's needs. The 
aerodynamic coefficients of the configuration are 
the least sensitive and converge first. A later 


convergence is obtained when the surface-pressure 
distribution is used as the convergence indicator. 
The most sensitive convergence test is the shape of 
the wake, which may continue to change long after 
the other parameters have converged. A typical 
convergence process of the aerodynamic coefficients 
of a cruciform missile configuration is shown in ■ 
Fig. 8 for the wings in a "plus" or in an x con- 
figuration. The convergence is fast. Only four 
iterations are needed for the aerodynamic coeffi- 
cients to converge, in spite of the large error In 
the results of the first iteration. The wake in 
these computations converged after seven iterations 
Agreement of the converged pitching-moment coeffi- 
cient with the data of R. Arieli (unpublished 
experimental results, Technion-IIT, Haifa, Israel, 
1982) is good, whereas the normal-force coefficient 
is within 5% of the measured value. This differ- 
ence will be discussed later. 

Since the changing shape of the wake deter- 
mines the solution of the nonlinear problem, it is 
safe to assume that the parameters of the wake and 
of its numerical treatment may influence the conver 
gence rate, or even determine if convergence is 
possible at all. These parameters may also influ- 
ence the converged results if the solution is not 
unique. 

The wake parameters that may be involved and 
have been studied here are 

1) The integration method of the free-vortex 
trajectories (Eq. (3)) 

2) The cutoff distance for the interaction of 
two close free vortices or for the interaction of a 
free vortex with a very close source panel 

3) The paneling method and its effect on the 
free-vortex trajectories 

4) The length of the integrated rolled-up wake 
(the distance to the effective far wake) 

5) An initial guess of the shedding angle of 
the free vortices 

6) The grid refinement in the wake and on the 

solid surfaces (1 

All three convergence criteria were considered. A 
convergence criterion of one count in the lift 
coefficient would suffice for engineering purposes. 
A more demanding convergence criterion would be a 
maximum change in the local pressure of 0.01 p m . 
Finally, changes between iterations in free- 
vortex trajectories should be limited to a maximum 
of one-tenth of the smallest vortex-panel dimension 
The criterion ACl < 0.0001 covered all of the 
above and was used as a single convergence 
criterion. 

Free-Vortex Trajectory Integration 

Equation (3) can be integrated by several 
methods. The general second-order integration of 
dy/dx = f(x,y) (the two-dimensional example was 
chosen for the sake of simplicity) would be 

y 1+ i = y t + R(Ax)|(l - io)f (x^y^ + ufj^ + f-j- . y ± 
+ f(x i' y i ) ]} + 0(Ax3) (8) 
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where R is an underrelaxation coefficient 
(0 < R < 1). The usual second-order integration 
uses u = 0.5. With id = 0, Eq. (8) is reduced to 
the simple, first-order Euler integration method, 
and u> = 1 gives an improved Euler scheme. In 
addition to these three methods, a second-order 
Runge-Kutta predictor-corrector was used. 

The simple Euler method was the least time 
consuming per iteration, whereas the Runge-Kutta 
method was, of course, the most time consuming. The 
differences between the converged solutions were too 
small to justify the increased computation costs 
(see also Refa. 8, 14, 20, 22). Therefore, the 
first-order Euler scheme is recommended. The con- 
vergence rate of this scheme in a scalar computation 
can be accelerated when the most recently updated 
wake shape is used to compute the induced velocities 
in Eq. (3) for the updating of the next grid point. 
This, however, prevents vectorization on a vector 
computer and the more efficient method has to be 
determined by trial and error. 

The Euler integration can be performed along a 
single vortex line out to the far wake, or can be 
marched downstream across the wake, correcting the 
first segments of all free vortices first, then the 
second segments, etc. No difference was found 
between the converged results of the two methods. 

An underrelaxation (R < 1) can be used to pre- 
vent the solution from diverging because of large 
fluctuations in the wake shape. Underrelaxation was 
usually not necessary in the present work, but 
Almosnino 25 reports that strong underrelaxation 
(0.3 < R < 0.5) had to be used for convergence of 
wakes shed from bodies. 

In conclusion, the integration method of the 
wake shape does not affect the uniqueness of the 
solution once it converges. 

Cutoff Distances for a Strong Interaction 

When a dense vortex lattice is used to improve 
the solution (discussed later), free vortices are 
apt to pass close to each other, at least during the 
iterations for the wake shape. The very high 
velocities induced at a close proximity of two 
potential vortices may cause the solution to diverge. 
This problem can be eliminated by preventing the 
induced velocity from increasing above a certain 
level. Three methods have been introduced in the 
past to achieve this purpose 14 : a "viscous core" 

(solid-body rotation) can be added to every vortex, 
two vortices can be combined into one when they are 
too close, or a constant induced velocity can be 
used instead. All three methods depend on an arbi- 
trary specification of a lower bound on the dis- 
tance separating two vortices. This distance is 
often called the "cutoff distance." Reference 20 
presents an interesting argument for the choice of 
a cutoff distance that depends on the parameters of 
the vortex model and on the free-stream velocity. 
However, it.'is based on empirical consideration and 
has no physical justification. 

In a wing-body interaction solved by the pres- 
ent method, 21 free vortices. that are shed from the 
lifting surfaces may pass too close to a bound vor- 
tex on another surface, or' to a source-panel on the 
body surface. As in the previously discussed case, 
the resulting strong interaction may have a cata- 
strophic influence on the iterative computation. A 


cutoff distance, similar to that above, must also 
be applied here. Without a sound theoretical basis 
for the specification of the cutoff distance, one 
has to rely on trial and error. Results of experi- 
mentation with the cutoff distance 22 indicate that 
it should be bounded from below by one-tenth of the 
smaller of the surface-panel dimensions, and from 
above by one-quarter of the same dimension. Above 
the upper bound the interaction is not fully 
accounted for. Below the lower bound the induced 
velocities are too strong and the solution may 
diverge. The wake roll-up in both cases is incor- 
rect. This means, however, that the solution is 
no longer unique and now depends on this cutoff 
distance. 

Effects of Paneling Scheme 

The paneling scheme (Fig. 6, trapezoidal, 
rectangular, and triangular) did not affect the 
results of the linear solutions (e.g., Fig. 7) 
which were in good agreement with experimental data 
obtained at' low angles of attack. 5 This was found 
not to be true in the nonlinear case. Whereas all 
the methods that accounted for leading-edge separa- 
tion could predict the lift coefficient fairly well 
(Fig. 9), the present method underpredicted the 
pitching moment of Ref. 26 (Fig. 10) because the 
pressure distribution (both in the chordwise and 
spanwise directions) was incorrect (Fig. 11). |The 
inboard pressures were too high and the outboard 
pressures were too low. The leading-edge suction 
peak was always too close to the leading edge when 
the Rom-Zorea paneling was used. 8 ’ 14 ’ 18 It is 
interesting to note here that, in spite of the 
locally incorrect pressures, any line integral of 
the pressures gave good results. A chordwise inte- 
gral resulted in a good prediction of the local 1 
section lift coefficient, and a spanwise integral 
gave a good approximation of the local contribution 
to the wing-root bending moment 22 (Fig. 12). The 
center of rotation and the center of gravity of the 
free vortices, as computed in the cross-flow plane 
(Fig. 13), did not agree with the experimental data 
of Ref. 27. 

These problems did not plague the delta-wing 
solutions of Ref. 15. A careful comparison of these 
solutions with the present ones showed that the 
vortices in the triangular panels along the leading 
edge (Fig. 1) in the present method were signifi- 
cantly weaker than those of Ref. 15. This could 
only be the result of the paneling method. 

In the present method the bound vortex lies 
completely within the triangular panel along its 
1/4-chord line, with the collocation point In .the 
middle of the 3/4-chord line (Fig. 14). This 
results in a bound vortex that is longer than those 
of the rectangular panels and results in a rela- 
tively shorter distance to the collocation point, 
both results contributing to a lower vortex inten- 
sity (Tj) in the solution of Eq. (4). These length 
differences become more pronounced as the leading- 
edge sweep angle is increased. In the paneling of 
Ref. 15 (Fig. 14), the bound vortex started at the 
inboard 1/4-chord of the triangular panel and' con- 
tinued out of the wing and into the flow, normal to 
the leading edge, ending on the chordwise line 
drawn through the outboard corner of the panel t 
(Fig. 14). Thus, it was of about the same length 
as the vortices in the rectangular panels with 
approximately the same distance to the; collocation 
point. The shorter vortex length resulted ima 
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higher vortex strength than in the present method 
and better pressure' distributions. 15 When the wing 
aspect ratio is increased and the leading-edge 
sweep angle is decreased, the paneling method of 
Ref. 15 increases the distance to the collocation 
point and overpredicts the strength of the leading- 
edge vortex. It seems that maintaining a constant 
ratio of the length' of the .bound vortex to its dis- 
tance to the collocation point would achieve the 
best results. This has been pointed out recently 
by Gordon and Rom. 28 In their calculations, the 
triangular panel was practically replaced by a 
standard rectangular panel and the computed results 
agreed well with experimental data. However, the 
physical meaning of- extending the leading-edge 
bound vortices into the outer flow field is not 
clear. In fact, the mathematical model of the 
delta wings of Refs, 15 and 28 resembles .a wing 
with sawteeth, which must certainly have different 
aerodynamic characteristics than delta wings . 

In conclusion, the results of the nonlinear 
mathematical model are not unique, since they depend 
on the vortex model and on its geometrical relation 
to the paneling scheme. 

The Length of the Rolled-up Wake 

The detailed .computation of the shape of the 
rolled-up wake (the Integration of Eq. 3) is per- 
formed from the vortex shedding points to a pre- 
scribed distance (S. w ) that is assumed to be the 
beginning of the "far wake." Further downstream the 
vortices are continued to infinity as straight 
lines. This practice saves computation time, since 
longer rolled-up wakes use more computer time and 
memory volume. 

- The effect of the prescribed distance to the 
far wake oh the computational results was small. 

The variation of the lift and .pitching-moment coef- 
ficients of. two delta wings (aspect ratios of 2.0 
and 1.0 at a = 20°) with (i w ) is shown in Fig. 15. 

A' rolled-up ^vake length .of half the mean aerodynamic 
chord was sufficient for the convergence of the 
aerodynamic coefficients. 

Initial Shedding Angle 

An initial guess of the wake shape has to be 
used in the first iteration cycle of the computation 
for the determination of the influence .coefficients 
of the free' vortices in Eq.,<5). Unless a specific 
wake shape from a previous solution is used, the 
simplest and easiest choice is of a linear wake 
(i.e., semi-infinite straight lines that leave the 
shedding points) deflected up or down from the free- 
stream direction at an angle (p). Angles between 
P = 0° (aligned with V„) and u = a (aligned with 
the chord) are used. A common choice 18 is p = a/2. 

Investigation of the influence of the initial 
shedding angle (p) value on the convergence charac- 
teristics of the method showed that p (within the 
above limits) had no effect .on the converged values 
of the aerodynamic coefficients (Fig. 16), and had 
but little effect on the convergence rate (Fig. 17a), 
even when the wake was initially .deflected .upward, 
too far away from its final position (Fig. 17b). 

Grid Refinement 

In a finite-difference solution of a differen- 
tial equation, a mesh refinement should, in princi- 


ple, converge to the solution of the differential 
equation itself, whereas in VLMs there is no proof 
that it should. The tangency boundary condition 
•is satisfied at the .collocation points only, while 
everywhere else the fluid "leaks" through the sur- 
face. This leak will not be eliminated by grid 
refinement. Still, the character .and trends of the 
variation of the solution with vortex-grid refine- 
ment is indicative of the uniqueness of the 
solution, 

- i 

The nonlinear VLH has two supposedly indepen- 
dent length scales. One is the length (Ax w ) of the 
free-vortex segments in the wake that is used in 
the integration of the wake shape (Eq. (3)). The 
other is the length (or chord) of the vortex panel 
(Axp), A reduction in (Ax„) should result in a 
smoother vortex line that is closer to the exact 
solution of the streamline equations. The effect 
of the reduction in the segment length on the lift 
coefficients of two delta wings (AR of 1.0 and 2.0 
at a = 20°) is shown in Fig. 18 for several con- 
stant vortex-panel sizes. When the length of the 
vortex panel AXp = c r /n c is kept constant, the 
ratio of panel length to vortex-segment length 
c r / (n c • Ax„) goes to infinity as AXy goes to zero. 

( 

The lift coefficient does exhibit a convergent 
character when the segment length is reduced. How- 
ever, when Ax w becomes much smaller than the panel 
length (typically less than one-quarter of Ax p ) 
the iteration process becomes unstable and diverges. 
Very long wake segments (typically more than five 
times longer than the panel) usually destabilize the 
computation. The differences between the solutions 
■with different segment lengths can be large 
(Fig, 18), The experimental lift coefficients for 
the same wings are shown for comparison in Fig. 18. 

Furthermore, Fig. 18 shows a very strong .effect 
■of the panel size on the lift coefficient. This 
effect is Investigated in Fig. 19, where the lift 
coefficients of the same two wings are plotted as a 
function of the total number of vertex panels on 
one-half of the wing (one-half was sufficient for a 
.symmetric solution). This is done for several con- 
stant wake segment lengths. The experimental lift 
coefficients are also shown. The lift coefficient 
.again shows a tendency to exponentially approach a 
constant value when the grid is refined. This 
final value depends, however, on the length of the 
free-vortex segment. 

To isolate the dependence .of the results on 
the two length scales, the previous results are 
replotted In Fig. 19 as the variations of the lift 
coefficient, the induced drag coefficient, and the 
location of the center of pressure with the increas- 
ing number of vortex panels for several constant 
ratios .of the free-vortex segment length to panel 
length. By -using this ratio as a parameter, the 
grid refinement applies simultaneously to both 
characteristic lengths. The lift and drag coeffi- 
cients show a tendency to reach their final values 
at or above 171 vortex panels (obtained by a divi- 
sion of the root-chord and semi-span into 18 equal 
parts). The final value of the lift coefficient 
depends on the wake-segment-to-panel-length ratio, 

A similar behavior .observed by Gordon and Rom, 26 
who used a different paneling scheme, indicates 
that this effect of the length ratio is not a 
result of the paneling method. 
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A comparison of the computed lift coefficients 
with the experimental data for three delta wings 
(Fig. 19) offers an approximate engineering formula 
for the choice of a length ratio that, when used 
for delta wings in the present method, results in 
convergence (resulting from grid refinement) to lift 
coefficients that are in fair agreement with the 
experimental data (for the aspect-ratio range tested 
here from AR = 0.5 to AR = 2.0). A ratio of 
(AXy) / (Ax p ) S (9 + AR 2 ) 1 / 2 /AR (or the simplified 
AXy/AXp s 3/AR) is recommended, although no reason- 
able explanation can be offered for this correlation. 

It is interesting to note that the correct 
location of the center of pressure (Fig. 14) is 
obtained even with a coarse vortex lattice and does 
not depend on the length-scale ratio. Apparently, 
the variation of the load distribution over the 
wing with both grid refinement and length-scale 
ratio is self-similar, so that the lift and pitching 
moment change simultaneously without moving the 
center of pressure. 

Another indication of convergence with grid 
refinement is the ratio (^/(Cl tan a), where Cp^ 

is the induced drag coefficient calculated by inte- 
gration of the forces that are exerted on the bound 
vortices by the induced cross-flow velocities. The 
coefficient Cp^ should, in principle, be equal to 

C L tan a and the ratio Cp^/^L tan a) should con- 
verge to 1.0. Grid refinement does indeed lead to 
convergence (Fig. 14), but only to a value of about 
0.9. Note that this ratio is absolutely independent 
of the length-scale ratio.. / 

Conclusions 

The convergence characteristics of the non- 
linear VLM have been thoroughly investigated. It is 
shown in this paper that the iterative solution 
procedure usually converges. Divergence can be pre- 
vented by a proper choice of a cutoff distance for 
the interaction of close singularities and by an 
underrelaxation of the iterative correction of the 
wake shape. :■ ' 

The computed results are independent of the 
specific integration method of the trajectories of 
the free vortices in the wake. The far-wake approx- 
imation can be applied at the relatively short dis- 
tance of half the mean aerodynamic chord behind the 
wing trailing edge with no detrimental effects on 
the computed results. The converged values of the 
aerodynamic coefficients are independent of the 
initial guess for the vortex shedding angle. 

Although grid refinement leads to convergence, 
the solution is not unique in this respect. The 
lift coefficient converges to different values, 
depending on the ratio of the free-vortex segment 
length to the panel length. An optimum length- 
scale ratio is recommended for engineering use. 

The local details of the solution, such as the 
load distribution, the shape of the wake, and the 
induced crossflow velocities have a strong depen- 
dence on the wing paneling method and vortex model. 
Thus, the solution is not unique in this respect. 

This nonuniqueness would nbt affect users who 
require only the aerodynamic coefficients of a 
single lifting surface. However, users who need 
the local details of the flow field or ones who use 


the method to compute interactions between several 
lifting surfaces and/or bodies, where the correct 
wake geometry is essental, should be aware of the 
limitations of the VLM and should ensure that the 
vortex model and paneling they use are adequate for 
the problem. 
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DATA: CONFIGURATION GEOMETRY 
FREE-STREAM VELOCITY, 

ANGLE OF ATTACK, a 
NUMBER OF SOURCE PANELS, N s 
NUMBER OF VORTEX PANELS, N v 
FREE-VORTEX SEGMENT LENGTH, Ax w 
LENGTH OF ROLLED-UP WAKE, S w 
INITIAL SHEDDING ANGLE, p 



Fig. 2 Flow chart of the iterative computation 
procedure. 
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Fig. 4 Pressure distribution on an ogive-cylinder 
at a = 6°. 
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Fig. 5 Relative error in the chordwise load dis- 
tribution on a two-dimensional flat-plate airfoil 
section, a = 4°. 


Fig. 3 The pressure distribution on an ellipsoid — 
differences between the numerical and analytical 
solutions . 
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Fig. 8 Typical solution — convergence process, a) Wings in a "plus" (+) position. 

b) Wings in an *V position. 
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Fig. 13 Computed cross-flow velocity field. 


14 





DELTA AR = 2.0 


P.deg 
A 0 
□ 10 
O 20 


M, deg 

O -20 OPEN SYMBOLS C NQR 

□ -15 SOLID SYMBOLS C M 

A -4 

25 T a = 20° (R. ARIELI, UNPUBLISHED 

( EXPERIMENTAL RESULTS, 1982) 


DELTA AR = 1.0 


25 20J 


® ~ C M C NOR 

15 * 15 


■ ■ ■ 


AR = 1.0 


Q □ 


Q AR = 20 

S g Q 

_J I I 

2 3 4 

NUMBER OF ITERATIONS 



1 ' MOMENTS REFERENCE POINT 

0 

... DIGITATED WINGS (CANARD +, WING x) 
IM 1 I I — : 

1 2 3 4 5 6 7 

NUMBER OF ITERATIONS 


I'ig. 17 Effects of the initial shedding angle on the convergence rate, a) Shedding angles between the 
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